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PACS 03.75.Lm - Tunneling, Josephson effect, Bose-Einstein condensates in periodic potentials 

PACS 64 . 70 . Tg - Quantum phase transitions 

PACS 67 . 85 . H j - Bose-Einstein condensates in optical potentials 

Abstract. - We investigate systems of bosonic particles at zero temperature in triangular and 
hexagonal optical lattice potentials in the framework of the Bose-Hubbard model. Employing 
the process-chain approach, we obtain accurate values for the boundaries between the Mott insu- 
lating phase and the superfluid phase. These results can serve as reference data for both other 
approximation schemes and upcoming experiments. Since arbitrary integer filling factors g are 
amenable to our technique, we are able to monitor the behavior of the critical hopping parameters 
with increasing filling. We also demonstrate that the g-dependence of these exact parameters is 
described almost perfectly by a scaling relation inferred from the mean-field approximation. 
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Introduction. — Over the last ten years ultracold 
atoms in optical lattices induced by standing waves of laser 
radiation have become an outstandingly important and 
Intensely studied testing ground for quantum many-body 
physics HJH]. Great prospects offered by these systems 
stem from the chance to investigate condensed-matter 
phenomena by simulating paradigmatic model Hamiltoni- 
ans in the laboratory [3]. In particular, the Bose-Hubbard 
Hamiltonian 4,5: has attracted a lot of attention, since it 
describes ultracold bosonic atoms in an optical lattice po- 
tential fairly well. This system exhibits a quantum phase 
transition from a superfluid to a Mott insulator upon in- 
creasing the lattice depth [6,7 . Its extensions even show 
further interesting phases, e.g. a supersolid state [8], when 
admitting particle-particle interactions between neighbor- 
ing sites [5] or introducing Bose- Fermi mixtures [TTJ] . 
' So far, most studies dealing with the Bose-Hubbard 
model have considered a square or a cubic lattice. For 
these particular lattice geometries the supcrnuid-insulator 
phase boundary has been calculated by various methods, 
such as mean- field approaches [4l [TlTH5] . the quantum ro- 
tor approach [16| , or a variational cluster formulation |17j . 
Arguably, the most precise results have been achieved by 
the strong coupling expansion [PEH2"0] and by Quantum 
Monte Carlo simulations [2TJ[22] for low filling factors of 
the lattice, and by means of the process-chain approach 
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for arbitrarily high integer filling |23j . 

Quite recently, the successful experimental realization 
of planar triangular and hexagonal lattices has been re- 
ported [21] ■ However, reliable theoretical data for the 
phase boundaries pertaining to these lattice types still 
seem to be missing, except for the single case of a triangu- 
lar lattice at unit filling (g = I), which has been covered 
by a strong coupling expansion [25j . Apart from the need 
to compare experimental results to accurate theoretical 
predictions, precise knowledge of the critical values of the 
hopping parameters would also be of great value to aid the 
development of new approximation schemes, and of future 
numerical methods. 

In this contribution we provide the phase diagrams 
for the Bose-Hubbard model with planar triangular and 
hexagonal lattice geometries. These two lattice types are 
depicted schematically in fig. Q] The process-chain ap- 
proach [26] in combination with the method of the effec- 
tive potential [27] [28] enables us to compute the phase 
boundaries with high precision, as has been demonstrated 
previously for square and cubic lattices [2"51l2l)] . 

For self-consistency, we start with a brief description 
of the Bose-Hubbard model, and give a short explana- 
tion of both the process-chain approach and the method 
of the effective potential, which provides the signature of 
the phase transition. We then present our results for the 
phase diagrams arising from triangular and hexagonal lat- 
tices, and state the corresponding critical values of the 
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Fig. 1: Schematic illustration of the triangular (a) and the 
hexagonal (b) lattice. Circles represent lattice sites occupied 
by bosonic particles; lines represent nearest-neighbor couplings 
due to tunneling processes between adjacent sites. 



hopping parameter (J/U) c and of the chemical potential 
(n/U) c . Since we can treat lattices with an arbitrary num- 
ber of particles per lattice site, i.e., with an arbitrarily 
high integer filling factor g, we are able to reveal that the 
critical values (J/U) c can be scaled such that they become 
(almost) independent of the filling factor. 

The model. — We study the homogeneous Bose- 
Hubbard model, given by the Hamiltonian 



H = — Y^ "*("» - !) ~ J Yl ^i ~ ^ ^ " l ' 
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which embodies in an elementary way the competition be- 
tween the kinetic energy due to tunneling processes and 
the potential energy associated with the repulsive interac- 
tion of bosons on the same lattice site. The operators c^ 
and a] are the bosonic annihilation and creation operators 
at site No. i, and n, is the corresponding number operator. 
We examine the case of zero temperature, which permits a 
single-band description, such that one only needs to con- 
sider the lowest Wannier state at each site. Moreover, an 
on-site approximation is made here, assuming that only 
particles sitting on the same lattice site interact with each 
other, each on-site pair contributing the amount U to the 
total interaction energy. Hopping processes of the bosons 
are restricted to adjacent sites; their strength is quanti- 
fied by the matrix element J. The subscript (i,j) at the 
kinetic-energy sum indicates that this summation only in- 
cludes pairs of neighboring sites. For the homogeneous 
systems studied here, the chemical potential (i is constant 
throughout the lattice. 

When expressing all energies in multiples of the on-site 
pair interaction energy U, we arrive at the dimensionless 
Hamiltonian 

H EB = ^^2fk(fk-l)-J/U^2 al&j - n/U^fii (2) 

containing two parameters, the hopping parameter J/U 
and the scaled chemical potential /J./U. 

The existence of a quantum phase transition from a 
Mott insulator to a supcrfluid [4;.6] in response to an in- 



crease of the hopping parameter is made plausible by in- 
specting the limiting cases: When J/U ^> 1 one has an 
almost ideal Bosc-Einstein condensate with all particles 
occupying the zero-quasimomentum Bloch state. This cor- 
responds to a supcrfluid with all particles dclocalized, and 
phase fluctuations being suppressed. The supcrfluid phase 
is characterized by long-range order and non-zero com- 
pressibility, d (n) /dfi 7^ 0. In the opposite limit J/U <C 1 
hopping is prohibited and all sites are decoupled from each 
other, so that the Hamiltonian ([2]) becomes diagonal in the 
occupation number basis. Minimizing the on-site energy, 
one finds that an integer number g = N/M occupies each 
site, with N denoting the total number of particles, and 
M the number of lattice sites. This phase is characterized 
by reduced density fluctuations and incompressibility, i.e. 
d (n) /dfi = 0. The ground state \m) for J/U = simply 
is a product state of Fock states with g particles on each 
site, 
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(3) 



where |0) is the particle- free vacuum. When starting in 
the Mott-insulating phase and increasing J/U from zero 
to higher values for a given, fixed chemical potential ///£/, 
there is a value (J/[/) p b at which the excitation gap van- 
ishes, marking the entrance into the supcrfluid regime. 

The Method. — In order to determine these values 
(J/U)pb of the hopping parameter at the phase bound- 
ary we make use of the method of the effective poten- 
tial [27H29] , which requires to add source and drain terms 
of constant strength r\ and rf to the Bosc-Hubbard Hamil- 
tonian ([2]): 

H W (.V,V1 = Hbh + ]T (t?*o 4 + v4) ■ ( 4 ) 

i 

The corresponding grand canonical free energy 

F(J/U,r,,r,*) = M (f (J/U) + £c 2ll (J/C/)|r,| 2 " j (5) 
with expansion coefficients 

CO 

c2n(j/u) = Y / 4:hj/ur (6) 

then is Legendre-transformed into an effective potential 
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r(j/v,i>,r) = f - -h>\ 2 + -im 4 +wn ■ (7) 

c 2 c| 

Odd orders of r\ vanish in the expansion ((SJ) of the free 
energy because of the [/(l)-symmetry of the augmented 
Hamiltonian (J3J). The expansion parameter |^| 2 of the 
effective potential (JT]) serves as the order parameter; it is 
given by 



u \ l dF 

/*/ N l dF 



i/V ' 



(8) 
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The Legendre pair 77 and ip* obeys the identity 



dT 

dip* 



-V, 



(9) 



the complex conjugate of this equation connects 77* and 
tjj. Now the original Bosc-Hubbard Hamiltonian © is re- 
covered from the augmented Hamiltonian (T5|) by setting 
77 = 77* = 0, which means that the relevant values of ip 
and ip* are those which render the effective potential T 
stationary. For low hopping strengths J/U, when the sys- 
tem is in its Mott-insulating phase, the coefficient C2 in 
the expansions ([5]) and ([7]) is negative, whereas C4 is posi- 
tive, leading to a minimum of T at \tp\ 2 — 0. The order 
parameter \tp\ adopts a non-zero value in the superfluid 
phase, signaling long-range order. The phase transition 
therefore takes place at that value of J/U for which I/C2 
vanishes, so that the minimum of the expression ((7|) starts 
to deviate from \ip\ 2 = 0. The upshot is that the phase 
boundary ( J/U) v \, equals the radius of convergence of the 
series (O for the coefficient C2. 

The coefficients a 2 °f that series are calculated within 
the process-chain approach, which is based on a diagram- 
matic evaluation 26]2§\ of Kato's perturbation series (30j . 
The Kato formula for the nth-order energy correction ex- 
perienced by a nondegenerate unperturbed state \m) in 
response to a perturbation V reads 



£l n) = tr 



J2 S ai VS a2 VS a3 . . . s a ™VS a ™ +1 

{a e } 



(10) 



Here the sum runs over all sets of n + 1 non-negative in- 
tegers OLi which obey the constraint Y) , at = n — 1. The 
linking operators S a are defined by 



S a = 






\m)(m\ 



P (0) 



(E£>-E?>) 



(0)n 



for a = 

for a > ' ( n ) 



where |z) denotes the unperturbed "intermediate" cigen- 
states, and E\ the corresponding unperturbed eigenval- 
ues. Kato's trace formula ([TU|) can be rewritten as a sum 
of matrix elements of the state \m) considered, 



(m\VS ai VS a2 . . . S^-'Vlm) 



(12) 



The number of such matrix elements (Kato-terms) quickly 
increases with the order n of perturbation theory In first 
order, the only Kato-term is (m\V\m), while n = 2 leads 
to (m\VS 1 V\m). These are precisely the well known first- 
and second- order energy corrections, as becomes obvious 
when inserting S 1 from eq. (jTTJ) . Each Kato-tcrm (|T2| can 
be viewed as a (sum of) closed process chain(s) consisting 
of n processes caused by the perturbation V, leading from 
the state \m) over various intermediate states \i) back to 
\m) again. When dealing with a homogeneous lattice sys- 
tem, many process chains can be combined into diagrams 
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Fig. 2: Diagrams of order v — 0, . . . , 3 in the hopping para- 
meter J/U for a triangular and for a hexagonal lattice. The 
symbol • indicates the creation of a particle, — ► corresponds to 
tunneling between two adjacent sites, and x to an annihilation 
process. Upper figure (a): Diagrams for the triangular lattice; 
weight factors are 0.) 1; 1.) 6; 2.) 6 and 30; 3.) 6, 138, 12, 30 
and 30. Lower figure (b): Diagrams for the hexagonal lattice; 
weight factors are 0.) 1; 1.) 3; 2.) 3 and 6; 3.) 3, 12, 6 and 6. 



by appending an appropriate weight factor. This proce- 
dure drastically reduces the numerical effort. A more de- 
tailed description of the application of this process-chain 
technique to the Bose-Hubbard model is given in ref. [55] . 
In our case, the unperturbed part of the Hamiltonian is 
site-diagonal, reading 



Ha = -= Y^ ««("i - 1) - l l / u E ; 



(13) 



The perturbation is given by the tunneling operators in 
combination with the source and drain terms artificially 
introduced in eq. (jl|): 



V 



j/u ^ ajoj- + Y yi* h i + va. 



(i,i) 



(14) 



Instead of using Kato's formulation for computing the to- 
tal energy corrections, we employ it for calculating the co- 
efficients a 2 °f the series (© for C2 only; the searched-for 
phase boundary (J//7) p b then is determined in a second 
step as the radius of convergence of this series. Because 
C2 is the coefficient of \tj\ 2 in the expansion of the free 
energy (J3J), it is associated with exactly one creation and 
one annihilation event of a particle. Hence, for calculat- 
ing its coefficients a 2 one nas to evaluate only diagrams 
containing one creation (symbolized by a dot: •) and one 
annihilation process ( x ) , together with v tunneling pro- 
cesses (— >)• 

Both the precise structure of the diagrams and their 
weight factors are determined by the geometry of the 
underlying lattice. Figure [5] lists the diagrams of order 
v = 0, . . . , 3 in the hopping parameter J/U for a triangu- 
lar and for a hexagonal lattice, together with their respec- 
tive weight factors. Each "hexagonal" diagram shown in 
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2 


3 


4 


5 


6 


7 


8 


9 


10 


Triangular 
Hexagonal 


1 
1 


1 
1 


2 
2 


5 

4 


14 
9 


41 
18 


129 
39 


416 
80 


1398 
180 


389 


1260 



Table 1: Number of diagrams to be evaluated when calculating 
the phase boundary of the Bose-Hubbard model for a triangu- 
lar and for a hexagonal lattice to ^th order in the hopping 
parameter J/U, corresponding to the order v + 2 of Kato's 
perturbation series. 



(b) is topologically equivalent to a "triangular" one in (a), 
but when taking three hopping processes into account a 
circular diagram turns up in the triangular case which has 
no hexagonal counterpart. In higher orders of the hop- 
ping parameter the number of "triangular" diagrams even 
becomes much larger than that of the "hexagonal" ones, 
as table [T] documents: The increase of the number of dia- 
grams with the number v of tunneling processes is much 
more pronounced in the triangular case. As another con- 
sequence of the geometric variation, the weight factors of 
corresponding diagrams generally differ for the two lattice 
types. 

The numerical value of a diagram is determined by go- 
ing through all permutations of its individual constituent 
processes; for each permutation one has to evaluate those 
Kato- terms which match it. The outcome then is mul- 
tiplied by the weight factor of the diagram in question. 
Finally the contributions of all diagrams occurring in a 
given order of perturbation theory are summed to yield 
the desired quantity a\ ■ For example, when consider- 
ing the hexagonal lattice with v — 3 tunneling processes, 
four diagrams depicted in fig. [5] (b) have to be dealt with. 
Each one of these leads to up to 5! = 120 different se- 
quences of processes which have to be matched with 3 
Kato-terms. Evidently the computational effort increases 
rapidly with the number v of tunneling processes taken 
into account: Both the number of Kato-terms and the 
number of diagrams proliferates quickly; in addition, the 
number of process permutations grows factorially with the 
order n = v + 2 of perturbation theory. 

Results. — For each preselected value of the chemical 
potential fJ,/U, the corresponding coefficients a^ of the 
series (JB]) for C2 show an almost geometric behavior, for 
both the triangular and the hexagonal lattice. As outlined 
above, the boundary (J/U) p b between the Mott insulating 
and the superfluid phase is given by the lowest J/U for 
which this scries diverges. Thus, for delineating the phase 
boundary we determine its radius of convergence by means 
of d'Alembert's ratio test [3"T] : 



(J/U) 
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The required extrapolation v 



linear fit of the ratios a 
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» is carried out by a 
over l/i/; the desired 
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value (J/U) p b then is the point of intersection with the 



Fig. 3: (Color online) Mott lobes for the triangular lattice (up- 
per panel) and for the hexagonal lattice (below) for various 
filling factors g. The dashed horizontal line \x/U = g — 0.5 
marks the axis of particle-hole symmetry which appears in the 
limit of large g. 



ordinate. This procedure also gives access to the rela- 
tive error of (J/U) p b'- Varying the set of coefficients a^ 
employed for the fit (e.g., taking only v = 4, . . . , 8) yields 
slightly different results; such fluctuations quantify the un- 
certainty of the final data. Here we employ the coefficients 
v = 2, . . . , 8, leading to an estimated relative error of less 
than 1% for the triangular case, and about 2% for the 
hexagonal one. 

The phase diagrams for the two lattice types are plot- 
ted in fig. [3] in the fi/U vs. J/U-pl&ne, for various filling 
factors g. The critical values (fj,/U) c and (J/U) c , i.e. the 
chemical potential and the hopping parameter at the tip of 
the respective Mott lobe, are listed in tabled Our result 
for the triangular lattice with unit filling compares favor- 
ably to the previous finding of Elstner and Monien [25] : 
These authors have stated (J/U) c = 0.037785, whereas we 
obtain (J/U) c = 0.03759; the deviation of about 0.5% is 
well within the estimated error margin. Because the coor- 
dination number z t = 6 of the triangular lattice is twice as 
large as that for the hexagonal one, z\ x = 3, the "triangu- 
lar" critical hopping strength at unit filling is substantially 
lower — by a factor of about 2.3 — than the "hexagonal" 
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Triangular 


Hexagonal 


9 


W(/)c 


(J/u) c 


0*/u)c 


(J/u) c 


1 


0.384 


3.759E-02 


0.360 


8.628E-02 


2 


1.432 


2.214E-02 


1.418 


5.075E-02 


3 


2.452 


1.574E-02 


2.442 


3.606E-02 


4 


3.463 


1.222E-02 


3.455 


2.799E-02 


5 


4.469 


9.984E-03 


4.463 


2.288E-02 


10 


9.484 


5.222E-03 


9.481 


1.196E-02 


20 


19.492 


2.674E-03 


19.490 


6.125E-03 


40 


39.496 


1.353E-03 


39.495 


3.100E-03 


50 


49.497 


1.085E-03 


49.496 


2.486E-03 


100 


99.498 


5.453E-04 


99.498 


1.249E-03 


1000 


999.500 


5.477E-05 


999.500 


1.255E-04 


10000 


9999.500 


5.480E-06 


9999.500 


1.255E-05 



Table 2: Critical values (/j,/U) c and (J/U) c for various filling 
factors g. For locating the tip of the respective Mott lobe, fi/U 
has been varied in steps of 0.001. Relative errors of (J/U) c 
axe less than 1% in the triangular case, and about 2% for the 
hexagonal lattice. 




10 100 1000 10000 

filling factor g 




1 10 100 1000 10000 

filling factor g 

Fig. 4: (Color online) Scaled critical values {J/U)l c according 
to eq. (|16[) for the triangular (upper plot) and for the hexagonal 
lattice (lower plot) vs. the filling factor g. Note the very fine 
scale of the ordinate. 



one. On the other hand, despite the fact that the coordina- 
tion number of the triangular lattice coincides with that of 
the simple three-dimensional (3D) cubic lattice, the cor- 
responding critical hopping strengths differ appreciably: 
The cubic lattice yields (J/U) c ~ 0.0341 for g = 1, see 
refs. |22H23) . amounting to a deviation of approximately 
9% from the triangular-lattice value. Inspecting the Mott 
lobes in fig. [3[ one also confirms that the critical chemical 
potential (/j,/U) c tends to g— 0.5 with increasing filling fac- 
tor g, as expected from the particle-hole symmetry which 
emerges in the large-<?- limit. 

Figure [3] also illustrates that the critical values (J/U) c 
decrease with increasing filling factor g. As we have shown 
previously [33] , in the cases of the 2D square and the 3D 
cubic lattices the g-dependence of the exact critical val- 
ues is quite well captured by the mean-field expression [4] 
for (J/U) c , even though the numerical agreement of the 
mean-field solution with the exact data is only moderate. 
Thus, the scaled critical values 



(j/uy c c = Vgia + T) 




l 



(J/U) c 



165(5 + 1) 

(16) 
are almost independent of g. Here we demonstrate that 
this finding also applies to the triangular and to the hexa- 
gonal lattice by plotting in fig. 0] the scaled data for both 
cases. As testified by the rather fine scale of the ordi- 
nate these scaled data are practically constant, with their 
residual variation amounting to only about 0.1%, which is 
an order of magnitude smaller than the estimated relative 
error committed in our present process-chain calculation. 
Finally, fig. [5] shows the triangular-lattice Mott lobes after 
applying the scaling (fTfj)) not only to (J/U) c , but to the 
entire phase boundaries. The scaled boundaries associated 
with different filling factors are quite similar; the remain- 



ing differences can be traced mainly to the particle-hole 
asymmetry of the Bosc-Hubbard Hamiltonian. Naturally, 
this asymmetry is reduced with increasing g. 

Conclusion. — We have presented fairly accurate 
phase boundaries for the homogeneous Bose-Hubbard 
model at zero temperature on both a triangular and on 
a hexagonal planar lattice, for filling factors ranging from 
unity to values so high that particle-hole symmetry is 
practically restored. The calculation has made use of the 
process-chain approach [26j , which already had proven its 
high fidelity for simple cubic lattices [53, 29 . Our numer- 
ical results can serve as benchmark data for other theo- 
retical approaches, and guide upcoming experiments with 
ultracold atoms in triangular and hexagonal optical lat- 
tices [21] ■ Furthermore, we have shown that the mean- field 
scaling ([TB)) of the critical values (J/U) c renders these data 
almost independent of the filling factor for both lattice 
types considered here. This ^-independence of the data 
scaled in this manner thus appears to be a general fea- 
ture of the Bose-Hubbard model, without being restricted 
to particular lattice geometries, while the lattice-specific 
scaled values (J/U)l° themselves may warrant further de- 
liberations. 
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Fig. 5: (Color online) Scaled phase diagrams for the triangular 
lattice with various filling factors g. The scaling provided by 
eq. (|16|l stretches the different Mott lobes such that their tips 
fall at the common value (J/U)l c . The remaining differences 
are caused by the particle-hole asymmetry of the Bose-Hubbard 
Hamiltonian. Again, the dashed horizontal line fi/U = g — 0.5 
marks the axis of symmetry which shows up for sufficiently 
large g. 
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